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ABSTRACT 

We place additional constraints on the three parameters of the dark matter halo 
merger rate function recently proposed by Parkinson, Cole & Helly by utilizing Smolu- 
chowski's coagulation equation, which must be obeyed by any binary merging process 
which conserves mass. We find that the constraints from Smoluchowski's equation are 
degenerate, limiting to a thin plane in the three dimensional parameter space. This 
constraint is consistent with those obtained from fitting to N-body measures of pro- 
genitor mass functions, and provides a better match to the evolution of the overall 
dark matter halo mass function, particularly for the most massive halos. We demon- 
strate that the proposed merger rate function does not permit an exact solution of 
Smoluchowski's equation and, therefore, the choice of parameters must reflect a com- 
promise between fitting various parts of the mass function. The techniques described 
herein are applicable to more general merger rate functions, which may permit a more 
accurate solution of Smoluchowski's equation. The current merger rate solutions are 
most probably sufficiently accurate for the vast majority of applications. 
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1 INTRODUCTION 

In current cosmological theory the mass density of the Uni- 
verse is dominated by dark matter. The most successful 
model of structure formation is that based upon the con- 
cept of cold dark matter (CDM). In the CDM hypothe- 
sis dark-matter particles interact only via the gravitational 
force. Since the initial distribution of density perturbations 
in these models has greatest power on small scales, the first 
objects to collapse and form dark-matter halos are of low 
mass. Larger objects form through the merging of these 
smaller sub-units. Consequently, the entire process of struc- 
ture formation is thought to proceed in a "bottom-up", hi- 
erarchical manner. 

Clearly then, the rate of dark-matter halo mergers is 
an absolutely crucial ingredient in models of galaxy and 
large-scale-structure formation, from sub-galactic scales to 
galactic and galaxy-cluster scales. The Press-Schechter (PS) 
formalism (Press & Schechter 1974) has long provided a 
simple, intuitive, and surprisingly accurate formula for the 
distribution of halo masses at a given redshift over a large 
range of mass scales and for a vast variety of initial power 
spectra. 

An elegant paper by Lacey & Cole (1993) — and similar 
work by Bond et al. (1991) and Bower (1991) — extended the 
work of Press and Schechter to determine the rate at which 



halos of a given mass merge with halos of some other mass. 
In addition to providing valuable physical insight, these 
merger rates have extraordinary practical value, having been 
applied to galaxy-formation models, e.g., if galaxy mor- 
phologies are determined by the merger history (Gottlober, 
Klypin & Kravtsov 1991); AGN activity (Wyithe & Loeb 
2003); models for Lyman-break galaxies (Kolatt et al. 1999); 
abundances of binary supermassive black holes (SMBHs) 
(Volonteri, Haardt & Madau 2002); rates for SMBH coales- 
cence (Milosavljevic & Merritt 2001) and the resulting LISA 
event rate (Menou, Haiman & Narayanan 2001; Haehnelt 
1994); the first stars (Santos, Bromm & Kamionkowski 2002; 
Scannapieco, Schneider & Ferrara 2003); galactic-halo sub- 
structure (Kamionkowski & Liddle 2000; Bullock, Kravtsov 
& Weinberg 2000; Benson et al. 2002; Somerville 2002; Stiff, 
Widrow & Frieman 2001); halo angular momenta (Vivit- 
ska et al. 2002) and concentrations (Wechsler et al. 2002); 
galaxy clustering (Percival et al. 2003) ; particle acceleration 
in clusters (Gabici & Blasi 2003); and formation-redshift 
distributions for galaxies and clusters and thus their distri- 
butions in size, temperature, luminosity, mass, and velocity 
(Verde et al. 2001; Verde, Haiman & Spergel 2002). 

However, as first noted by Lacey & Cole (1993; see 
their Section 3.1) and as we demonstrated in Benson, 
Kamionkowski & Hassani (2005; hereafter Paper I]) these 
merger rates are fundamentally flawed. The extended-Press- 
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Schechter (ePS) formulae for merger rates are mathemati- 
cally self-inconsistent, providing two different results for the 
same merger rate, as was also pointed out by Santos, Bromm 
& Kamionkowski (2002), which are increasingly discrepant 
for larger mass ratios. This ambiguity will be particularly 
important for, e.g., understanding galactic substructure and 
for SMBH-merger rates (Erickcek, Benson & Kamionkowski 
2006). Even the smaller numerical inconsistency for merg- 
ers of nearly equal mass may be exponentially enhanced 
during repeated application of the formula while construct- 
ing merger trees to high redshift. Moreover, the ambiguity 
calls into question the entire formalism, even when the two 
possibilities seem to give similar answers quantitatively*. 
Recently, Neistein & Dekel (2008b) have described an al- 
ternative derivation of merger rates from ePS theory and 
an implementation which appears to be self-consistent. Our 
goal here, however, is not to find a merger rate kernel which 
reproduces the results of ePS since that theory does not re- 
produce results found in N-body simulations. Instead, we 
aim to find a merger kernel which agrees well with N-body 
results and correctly evolves mass functions designed to fit 
the halo mass function measured from N-body simulations. 

In Paper I, we discussed the mathematical requirements 
of a self-consistent theory of halo mergers. As recognized al- 
ready (Silk & White 1978; Cavaliere, Colafrancesco & Menci 
1992; Sheth 1995; Sheth & Pitman 1997), the merger pro- 
cess is described by the Smoluchowski coagulation equation. 
Indeed, several authors have employed the Smoluchowski 
equation in precisely this way (Cavaliere, Colafrancesco & 
Menci 1991; Cavaliere, Colafrancesco & Menci 1992; Cava- 
liere & Menci 1993; Menci et al. 2002). This equation sim- 
ply says that the rate at which the abundance of halos of a 
given mass changes is determined by the difference between 
the rate for creation of such halos by mergers of lower-mass 
halos and the rate for destruction of such halos by mergers 
with other halos. The correct expression for the merger rate 
must be one that yields the correct rate of evolution of the 
halo abundance when inserted into the coagulation equa- 
tion. The problem is thus to find a merger rate, or "kernel," 
that is consistent with the evolution of the halo abundance, 
either the Press-Schechter abundance or one of its more re- 
cent N-body-inspired variants (Sheth, Mo & Tormen 2001; 
Jenkins et al. 2001). 

The apparent simplicity of the mathematical problem 
is in fact quite deceptive. The Smoluchowski coagulation 
equation is in fact an infinite set of coupled nonlinear dif- 
ferential equations. The equation appears in a variety of ar- 
eas of science — e.g., aerosol physics, phase separation in liq- 
uid mixtures, polymerization, star-formation theory (Allen 
& Bastien 1995; Silk & Takahashi 1979), planetesimals 
(Wetherill 1990; Malyshkin & Goodman 2001; Lee 2000), 
chemical engineering, biology, and population genetics — so 
there is a vast but untidy literature on the subject (although 
see Leyvraz 2003 for an illuminating review). It has been 

* Extended Press-Schechter theory discusses the trajectories of 
points in the primordial density field as the smoothing scale is 
reduced. It is the association of such trajectories with halo masses, 
which is not necessarily well-defined (see, for example, Porciaini, 
Dekel & Hoffmann (2002) who find that only 40% of proto-halo 
regions contain peaks of the density field), that leads to these 
problems with the derived merger rates. 



studied a little by pure and applied mathematicians (Aldous 
1999). Still, solutions to the coagulation equation are poorly 
understood. Furthermore, there is virtually no literature on 
the problem we face: i.e., how to find a merger kernel that, 
when inserted into the coagulation equation, yields the de- 
sired halo mass distribution and its evolution as a solution. 

In Paper I we solved Smoluchowski's equation numer- 
ically subject to two additional physically-motivated con- 
straints (or regularization conditions). First, the merger ker- 
nel should be positive for all masses^ and, second, the merger 
kernel should be a smooth function of its arguments. The 
second regularization condition was imposed by minimiz- 
ing the second derivatives of the kernel with respect to the 
two arguments. This allowed us to find unique solutions to 
Smoluchowski's equation for several power-law power spec- 
tra. These solutions were in reasonable agreement with the 
limited N-body data available for comparison. 

We have attempted to extend the work carried out in 
Paper I to the physically interesting case of a CDM power 
spectrum, using a much improved calculation able to span 
a wider dynamic range of halo masses. While we have been 
able to obtain solutions to Smoluchowski's equation for such 
power spectra we find that they do not correspond to the 
merger rates that occur in N-body simulations of structure 
formation, which we consider to provide the "correct" so- 
lution for our purposes. Specifically, while our solutions are 
correct solutions of the Smoluchowski equation, they do not 
reproduce the progenitor mass functions of halos as mea- 
sured from N-body simulations. (In fact, they typically per- 
form worse than the standard extended Press-Schechter the- 
ory.) This is, perhaps, not surprising as our results are con- 
trolled by the regularization conditions applied. We required 
solutions to be smooth, which seems reasonable, but our 
minimization is forced to pick the smoothest solution in or- 
der to find a unique answer. It seems, therefore, that the 
Universe uses a smooth merger rate, but not the smoothest. 
Trials with alternative regularization conditions have failed 
to find any more successful approach. 

Until such time as improved physical understanding of 
the physics governing the merger kernel is uncovered we have 
adopted a different approach to finding a suitable merger 
kernel. Recently, several authors have used the Millennium 
Simulation to provide constraints on either the merger rate 
(Fakhouri & Ma 2007) or the progenitor mass functions of 
dark matter halos (Neistein & Dekel 2008a). In particular, 
Parkinson, Cole & Helly (2008; hereafter PCH08) utilized 
progenitor mass functions to constrain the parameters of 
a function used to modify the standard extended Press- 
Schechter merger rate function and thereby constructed a 
merger tree binary split algorithm which accurately matched 
those progenitor mass functions. This approach is success- 
ful, but is limited by the accuracy and extent of the N-body 
data. In this paper, we demonstrate how we can use Smolu- 
chowski's equation to provide an additional constraint on 



t A negative merger rate could be considered to describe a spon- 
taneous fission process, but such processes should have a different 
dependence on halo abundances and so require the inclusion of 
additional terms in Smoluchowski's equation. We ignore such pro- 
cesses in this work. 
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such algorithms. This approach is powerful for the following 
reasons: 

(i) The PCH08 merger tree algorithm is a binary split 
algorithm and so should obey Smoluchowski's equation; 

(ii) Constraining the three free parameters of the PCH08 
algorithm is much easier than attempting to numerically 
invert Smoluchowski's equation subject to complicated reg- 
ularization conditions; 

(iii) Smoluchowski's equation applies over all halo masses, 
while even the Millennium Simulation has a limited dynam- 
ical range of masses that it can probe. 

The disadvantage of this approach is that it assumes a par- 
ticular functional form for the merger kernel (with free pa- 
rameters that are to be fit). There is no guarantee that this 
functional form will permit a solution to Smoluchowski's 
equation for any values of its free parameters and, in fact, 
we will show that it does not. Nevertheless, we can find the 
best-fit to the Smoluchowski equation. Further free param- 
eters could be introduced as constraints, of course, which 
should result in improved agreement with Smoluchowski's 
equation. 

The remainder of this paper is arranged as follows. In 
§2 we describe the Smoluchowski equation as used in this 
work and, in particular, how it is applied to the construction 
of merger trees. In §3 we present the results of fitting the 
parameters of the PCH08 model to N-body progenitor mass 
functions and to the Smoluchowski equation and examine 
how these influence the evolution of the dark-matter halo 
mass function. Finally, in §4 we examine the implications of 
these results. 



2 SMOLUCHOWSKI'S EQUATION AND 
MERGER TREES 



t The equation applies to any additive quantity conserved 
through the coagulation process. 



kM\M2- In this work, the objects we consider are dark mat- 
ter halos, and the mass is the total mass of those halos. The 
second analytic solution is of particular interest as it corre- 
sponds to the solution for dark matter halos in a Universe 
with a P(k) = k n power spectrum when n — 0. The other 
two analytic solutions are not cosmologically interesting. 

It is convenient to work with scaled variables for mass 
and time. As shown in Paper I, the natural time variable is 
t — — \n[S c )(t)] where 5 c (t) is the extrapolated linear the- 
ory overdensity required for halo collapse in the spherical 
top-hat collapse model (as appears in the Press-Schechter 
expression for the halo mass function). We also choose to 
express masses in units of the characterstic mass scale M* 
(defined such that a(M„[t]) = 5 c (t) = cxp(-r) where a 2 (M) 
is the fractional mass variance of the linear density field ex- 
trapolated to Z = (we use Z to indicate redshift so as to 
distinguish from the mass variable z). With these choices, 
the Smoluchowski equation becomes: 

r*/2 

y(z;r) = / q(z — z , z; r)n(z — z \ r)n(z'; r)dz' 
Jo 

f-OO 

- \ q(z,z'-T)n(z;T)n(z';T)dz', (2) 
Jo 

where z — M/M*(t), n(z) is the distribution of halo masses, 
y(z) = dn(z)/dT and q{z\, 22; t) is the merger rate function 
in these new variables. For the specific case of power-law 
power-spectra (which are, of course, scale-free) there is no 
explicit time-dependence with this choice of units and we 
can write 

j-z/2 

y(z) — / q(z — z' , z')n(z — z')n(z')dz 
Jo 

/>oc 

— / q(z, z')n(z)n(z')dz' . (3) 
Jo 

Thus, a single solution, valid at all times can be found for 
power-law power-spectra. For CDM power-spectra, which 
are not scale free, the merger kernel can in principle depend 
explicitly on time (although the choice of scaled variables 
will minimize this dependence). In principle, therefore, we 
could use Smoluchowski's equation at each point in time (r) 
to provide constraints on the merger rate function. We re- 
tain the variable choices z and r as we expect the solutions 
q(zi, 22; t) to be only slowly changing with r in these vari- 
ables, since the CDM power-spectrum is close to power-law 
over a wide range of masses. In practice, we will use only the 
r = ro (i.e. present day) Smoluchowski equation, although 
our methods could be easily applied to other redshifts if re- 
quired. 

It is well known that the Press-Schechter mass function 
is not a good description of the mass functions found in N- 
body simulations of CDM structure formation (Sheth, Mo 
& Tormen 2001; Jenkins et al. 2001). Our methods are ap- 
plicable to any mass function n(z) and so we can replace the 
Press-Schechter mass function with a fitting formula such as 
that from Sheth, Mo & Tormen (2001): 

nsMT(z;r) = \ -Asmtcx(z) - <yZ, T ^ > ex.p[-x' 2 (z, r)/2] 
\ n z z 

x(l + l/x' 29SMT ), (4) 
where x'(z,t) = ^/osmtx(z,t), x(z,t) = exp(-r)/a(z), 



In this section we will describe the techniques used in this 
work and the specific implementations used to construct pro- 
genitor mass functions and evolving dark-matter halo mass 
functions. 

Smoluchowski's equation describes the changing distri- 
bution of "masses" ^ of objects growing through coagulation. 
Given a distribution, n(M; t), of object masses M at time t, 
the Smoluchowski equation gives the rate of change of this 
distribution: 

h(M;t) = / Q(M - M',M';t)n{M - M';t)n(M';t)dM' 
Jo 

poo 

- / Q(M,M';t)n(M;t)n(M';t)dM', (1) 
Jo 

where Q(M\,M2\t) encodes the merger rate between ob- 
jects of mass Mi and M2 at time t. It is, therefore, sym- 
metric in its arguments, i.e. Q(M\,M2\t) = Q(M2, Mi;t). 
The first term in eqn. (1) represents creation events, while 
the second represents destruction events. Only three forms 
for Q(M\,M2\t) are known to permit analytic solutions of 
Smoluchowski's equation: Q(Mi, M2; t) = k (where A: is a 
constant), Q(M 1 ,M 2 ]t) = k(M 1 + M 2 ) and Q(M 1 ,M 2 ]t) = 
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flSMT = 0.707, gsMT = 0.3 and ^4smt(~ 0.3222) is chosen 
such that the mass density in halos equals the mean density 
of the Universe. 

Smoluchowski's equation involves the time derivative of 
this function, which is given by 



2/smt(z,t) = ii S mt(z,t) 



x' 2 (z,t) + 



2gsMT 



1 + x' 2 9SMT ( z , t) 



- 1 



.(5) 



As shown in Paper I, the standard extended Press- 
Schechter theory predicts a merger kernel 



Q e ps(Mi,M 2 ;t) 
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exp 
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(6) 



where we have adopted the notation 02 = o(M 2 ), etc. This 
is, of course, asymmetric in the two mass arguments. In 
scaled variables, this becomes 



q e ps(zi,z 2 ;T) = 



z\ cr 2 



<7f Zf 



din Of 



d In zt 



dlna 2 



din Z2 



(1 



exp 



*f>I. 
„-2r 



2\3/2 



(7) 



where, as shown in Paper I, the mean cosmic density, po, is 
removed through an appropriate choice of units. 

When computing the mass variance, o(M), we use the 
same power spectrum as used to generate the initial con- 
ditions for the Millennium Simulation since we will fit to 
conditional mass functions measured from that simulation. 
This power spectrum is integrated under a spherical top-hat 
real-space window function to obtain o{M). 



2.1 Application to Merger Trees 

A primary goal of this study is to provide accurate merger 
rates for use in the construction of dark matter halo merger 
trees. "Accurate" here is defined to mean that the merger 
trees so constructed should produce the correct distribution 
of progenitor halo masses at earlier times and, consequently, 
produce the correct evolution of the halo mass function. Ide- 
ally, this evolution should remain precise even when trees are 
constructed over large fractions of cosmic history. The need 
for such accurate trees is highlighted by the work of Benson 
et al. (2006) (see their Fig. 9), who were forced to restart 
their calculations of reionization at periodic intervals in red- 
shift due to accumulating inaccuracies in their merger trees. 
Constructing merger trees also provides a means of com- 
paring our results for the merger rate with measurements 
from N-body simulations. In this Section, therefore, we will 
describe how we construct merger trees. 



2.1.1 Merger Tree Construction Algorithm 

We adopt the algorithm described by Cole et al. (2000) and 
PCH08 to construct merger trees. This algorithm considers 



binary splits of the merger tree with a correction for accre- 
tion of mass in halos below the mass resolution of the tree 
and uses adaptive timesteps to ensure that the probability 
of non-binary mergers remains small. As Cole et al. (2000) 
note, their algorithm performs well due to the subtle choice 
of drawing a progenitor mass from the lower half of the ePS 
progenitor mass distribution function (i.e. for a parent halo 
of mass Mo they use only the < Mo/2 region of the pro- 
genitor mass function when computing the probability of 
a binary split). It is interesting to realize that this choice 
effectively symmetrizes the Smoluchowski merger kernel as 



?ePS,sym(-3l, Z 2 ) = 



9ePs(2l,22) 
9ePs(Z2,Zl) 



if zi < {zi + z 2 )/2 
if zi > (zi + z 2 )/2. 



(8) 



Other algorithms have been proposed for constructing 
merger trees (Kauffmann & White 1993; Cole et al. 1994; 
Somerville & Kolatt 1999)-we select the Cole et al. (2000) 
because it allows only binary mergers and is therefore consis- 
tent with a treatment via the Smoluchowski equation. We 
refer the reader to Cole et al. (2000) and PCH08 (their 
Appendix A) for a full description of the merger tree con- 
struction algorithm, but address the aspects of the imple- 
mentation unique to the present work below. 

Cole et al. (2000) compute two quantities, F(z,t) and 
P(z, t), which give the fraction of mass gained through ac- 
cretion and the probability of a binary split respectively for 
a halo of mass z at time r and during some time period 
St. In terms of quantities used in this work, the parameters 
F(z,t) and P(z,t) are given by: 



F{z,t) = St 



and 



z'R(z'; r)dz' , 



72 



P(z,r) = St / R{z';T)dz, 



where 



R(z') 



n(z'; r)n(z — z'; r)q(z' , z — z'; r) 
n(z;r) 



(9) 



(10) 



(11) 



Zmin = M m in/M* (t) and M m in is the mass resolution of the 
merger tree. We have expressed this rate in terms of the 
merger kernel in Smoluchowski's equation. PCH08 give ex- 
pressions for F and P in terms of the ePS expression for 
the mean number of progenitors of mass Mi expected in a 
small timestep. The two ways of expressing these functions 
are equivalent — we use the version depending on the merger 
kernel to make clear the connection to Smoluchowski's equa- 
tion as used in this work. We adaptively choose timesteps 
St to ensure that F <Si 1 and P <JC 1, to keep the possibility 
of multiple fragmentation small and to ensure that the halo 
mass can change by only a small amount due to accretion 
in any given timestep. 

Using the Press-Schechter mass function for n(z; r) and 
the ePS merger rate q c ps(zi, z 2 ;t) (cqn. 7) recovers the ex- 
pressions given by Cole et al. (2000) for F and P. Equa- 
tions (9) and (10) are more general however, allowing any 
mass function n(z; r) and merger rate function q(zi,z 2 ;T) 
to be used. In particular, PCH08, propose that R(z';t) be 
modified by multiplying by a function 



G(<7i/<7f,* c (T)/<7f) = Gofax/arP [* c (r)/ff f p 



(12) 
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where Go, 71 and 72 are free parameters. In terms of the 
Smoluchowski equation we can similarly multiply the ePS 
merger kernel by the same function G{ai/ai,8 c {T)/ai). We 
will call this modified merger kernel q' cPS and label such ker- 
nels as (Go, 71, 72)- However, we must also account for the 
fact that the PCH08 algorithm implicitly utilize the Press- 
Schechter mass function, while we actually want to repro- 
duce the evolution of the Sheth, Mo & Tormen (2001) mass 
function and must therefore use that when solving Smolu- 
chowski's equation. Appendix A describes how to correctly 
implement the Smoluchowski kernel in this case. 



random sampling. We impose a minimum and maximum 
number of halos to simulate from each decade of halo mass 
to ensure good sampling of the entire mass range of interest, 
and span a large enough range of initial masses to ensure 
that our results are unaffected by the necessarily limited 
range of masses probed. 

We then apply the Cole et al. (2000) method to evolve 
these halos backwards in time until a time T2 is reached. 
Summing the progenitor halo masses of all of the T\ halos 
allows us to construct the halo mass function at T2. This is 
then compared to the expected mass function hsmt(z; T2). 



2.1.2 Building Progenitor Mass Functions 

Utilizing the above algorithm we can, given a halo of mass 
Z{ at the present day, construct a merger tree back to some 
earlier time. By constructing a large number of such trees 
and averaging over their progenitors we can construct an es- 
timate of the progenitor mass function at that earlier time. 
Specifically, we wish to construct progenitor mass function 
to compare to those determined by Cole et al. (2008) from 
the Millennium Simulation. Cole et al. (2008) estimated pro- 
genitor mass functions of Z = halos by locating all Z = 
halos with a mass within a factor \f2 of some mass Mi and 
then finding all progenitors of those halos at redshifts 0.5, 
1, 2 and 4. To compare to these N-body data-* we first pro- 
duce a sample of halo masses M between Mf/^/2 and \f2M{ 
drawn from the Sheth, Mo & Tormen (2001) mass function. 
We then construct a merger tree for each such halo back to 
z — 4 and cumulate the progenitor halo masses to estimate 
a mean conditional mass function. 

We construct enough trees that our comparison is lim- 
ited by noise in the N-body results and not by the limited 
number of trees constructed. We then determine a \ 2 statis- 
tic using 



X 



/££(Mi|M f )-/£g(Mi,M f ) 



<f f (Mi|M f ) 



(13) 



where / cm f (Mi \M{ ) is the conditional mass function of halos 
of mass Mi which are progenitors of a halo of mass Mf at 
Z = 0, and superscripts refer to the Millennium Simulation 
(MS) and Monte-Carlo (MC) merger trees constructed using 
our algorithm. The values of cr^f f (Mi|Mf) used are deter- 
mined from the number of progenitor halos in the Millen- 
nium Simulation assuming Poisson statistics. This process 
is repeated for different values of the parameters G , 71 and 
72- 



2.1.3 Evolving the Halo Mass Function 

In order to study the evolution of the halo mass function, 
usmt(z;t), we begin by drawing a sample of dark matter 
halo masses at random from that mass function at some 
initial time n corresponding to Z = 0. We use a quasi- 
random method (specifically a Sobol sequence; Press et al. 
2007 section 7.7) to produce a non-uniform sample of masses 
while minimizing fluctuations in the mass function due to 



§ Available from 
http: //star-www. dur.ac .uk/~cole/merger_trees/MS_data/. 



2.1.4 Constraints from Smoluchowski' s Equation 

Given a merger kernel, QcPSj we can > f° r an y combination of 
parameters (Go, 71, 72), evaluate the creation and destruc- 
tion integrals in Smoluchowski's equation and thereby de- 
termine y(z) for that merger kernel. In general, a merger 
kernel will not provide a precise solution to Smoluchowski's 
equation. We can measure the ability of any given merger 
kernel to solve Smoluchowski's equation by evaluating the 
mean absolute difference between the resulting y(z) and that 
found by direct differentiation of the mass function with re- 
spect to time (which is, of course, the required solution to 
give the correct evolution of the mass function). Therefore, 
we compute 



iV 

«4e 



y(zi) - ysMT(zi) 



USMT(Zi) 



(14) 



where |/smt(z) is the rate of change of the Sheth, Mo &: 
Tormen (2001) mass function, and we average the absolute 
error over N — 16 different masses, z;, equally spaced in log z 
between z = 10 -5 and z = 10 3 . We have chosen to weight by 
1/tismt{z) so that we are comparing the rate of change of 
the mass function per halo. The choice of how to judge any 
merger kernel's degree of success in solving Smoluchowski's 
equation is somewhat subjective — the above choice ensures 
that we require the solution to be most accurate where the 
timescale for change in the mass function is most rapid (i.e. 
for the most massive halos) . 



3 RESULTS 

We find 
best-fit parameters of (Go, 71, 72) = (0.605,0.375,-0.115) 
and (0.775, 0.275, -0.225) for fits to the Millennium Simu- 
lation conditional mass functions and Smoluchowski's equa- 
tion respectively. Our best-fit parameters for matching con- 
ditional mass functions differ slightly from those obtained by 
PCH08 due to our choice to weight our goodness-of-fit mea- 
sure by the inverse Poisson errors on the N-body conditional 
mass functions. As expected, both constraints favour Go < 1 
(which reduces the overall merger rate) and 71 > (which 
boosts the ratio of low mass to high mass progenitors) — 
PCH08 discuss the reasons for these expectations. Figure 1 
shows constraints in the (Go, 71) plane (i.e. on slices of fixed 
72 in the 3-D parameter space), for both fitting to condi- 
tional mass functions (blue contours) and Smoluchowski's 
equation (red contours). In the left and right-hand panels 
we show results corresponding to the best-fitting values of 
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Figure 1. Constraints in the Go-71 parameter plane for two values of 72 (72 = —0.125 in the left-hand panel, corresponding to the 
best fit to Millennium Simulation conditional mass functions; 72 = —0.225 in the right-hand panel corresponding to the best solution of 
Smoluchowski's equation). The solid blue dot shows the best-fit values obtained by fitting to progenitor mass functions measured from 
the Millennium Simulation, while blue contours show the constraints obtained from this dataset. The red point shows the best-fit values 
obtained by fitting to the Smoluchowski equation, while the red contours map locii of constant mean error between y(z) and the fitted 
value as defined in eqn. (14). The constraints shown correspond to a slice through the 3-D parameter space at constant 72. The open blue 
circle shows the constraint obtained by PCH08 using the same Millennium Simulation data set. This differs slightly from our equivalent 
constraint due to our choice to utilize the measurement errors in the progenitor mass functions in our fit. Nevertheless, the two are quite 
similar as expected. (Note that PCH08 found a best fitting value of 72 = —0.01.) 



72 for conditional mass functions and Smoluchowski's equa- 
tion respectively. 

The conditional mass functions provide a good con- 
straint on (Go, 71) with a degree of degeneracy between 
the parameters. Smoluchowski's equation shows a much 
larger degeneracy between these two parameters and, in fact, 
favours regions which lie in a plane in the full 3-D parameter 
space. This plane is, however, thin, thereby strongly ruling 
out large regions of the parameter space. Remarkably, the 
degeneracy seen in fitting Smoluchowski's equation is similar 
to that from fitting conditional mass functions. More impor- 
tantly, the constraints from both methods intersect, permit- 
ting solutions which both match the conditional mass func- 
tions and are reasonable solutions to Smoluchowski's equa- 
tion. This is an important point, as it implies that N-body 
merger trees are consistent with a binary merger hypothesis. 

We could, in principle, combine these two constraints 
to obtain an overall best-fitting model. Such a procedure is, 
however, somewhat arbitrary. Firstly, while we have a statis- 
tically meaningful (in the sense that we can use it to assign 
relative probabilities to models) measure of goodness-of-fit 
for the fits to conditional mass functions, our goodness-of- 
fit to Smoluchowski's equation is not statistically meaning- 
ful. Any combination of the two constraints therefore re- 
quires some amount of judgment as to their relative im- 
portance. Furthermore, the (Go, 71, 72) modification to ex- 
tended Press-Schechter merger rates does not give a "good- 
fit" to conditional mass functions (it is certainly a much bet- 
ter fit than unmodified extended Press-Schechter, but clear 



systematic differences from the N-body data remain) and 
does not precisely solve Smoluchowski's equation. 

It is, however, obvious that any choice of a "best-fit" 
model to both constraints must lie somewhere in the plane 
defined by Smoluchowski's equation between the red and 
blue dots in Fig. 1. We will illustrate below the differences 
in results when these two best-fit values are obtained — any 
combined best-fit should show intermediate behaviour. 

It is instructive to examine how well our various mod- 
els solve Smoluchowski's equation. Figure 2 shows the dif- 
ference between the rate of change of the mass function 
per halo as a function of halo mass determined from our 
merger kernels and the value required for a correct solution 
of Smoluchowski's equation. The standard extended Press- 
Schechter merger kernel (green line) fails significantly over 
the entire range of masses plotted (particularly so at the 
massive end). The blue lines indicate results using merger 
kernels constrained to match the conditional mass func- 
tions from the Millennium Simulation. Over the range where 
these conditional mass functions provide a good measure of 
the merger kernel (approximately 0.001 < z < 10) the fit- 
ted merger kernels provide a significantly better solution to 
Smoluchowski's equation than does the unmodified extended 
Press-Schechter kernel. Note also that our fit (solid blue line) 
and that of PCH08 (dotted blue line) are very similar over 
this range. At z ;>10 the Millennium Simulation conditional 
mass functions provide only weak constraints on the merger 
kernel and, consequently, the kernels fitted to these func- 
tions provide worse solutions to Smoluchowski's equation 
in this regime (although still better than the original ex- 
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Figure 2. The rate of change of the halo mass function per 
halo, y(z)/n(z), relative to that obtained by direct differentia- 
tion of the Sheth, Mo & Tormen (2001) mass function (the hor- 
izontal black line indicates the true rate of change expressed in 
this way). Coloured lines show this quantity computed using the 
Smoluchowski equation using various merger kernels, i.e. using 
l'ePS wn "h different parameters (Go , 71 , 72)- The green line shows 
the result of using the standard extended Press-Schechter kernel, 
(1, 0, 0), while the red line shows the results for the merger kernel 
which best matches the true rate of change per halo as judged 
by eqn. (14). The solid blue line shows results for the merger 
kernel that best fits progenitor mass functions in the Millennium 
Simulation. The dotted blue line shows the same using the fit pa- 
rameters of PCH08. For all calculations, we use the Sheth, Mo &: 
Tormen (2001) mass function in Smoluchowski's equation. 



cates progenitor mass functions obtained using the merger 
kernel which best solves Smoluchowski's equation. This is 
clearly intermediate in success between an unmodified ex- 
tended Press-Schechter kernel and kernels constrained to 
match the conditional mass functions, as may be expected. 
(There are cases where it performs significantly better how- 
ever, for example, for the Z = 4 conditional mass function 
of W ZA hr x MQ halos.) 

Finally, in Fig. 4 we show total (i.e. not conditional) 
mass functions of halos from z — to z — 4. The left 
hand panel shows this function in physical units, while the 
right hand panel shows the fraction of mass in halos with 
v — S c (z) /a(M) per logarithmic interval of v (c.f. Figure 4 
of PCH08. In this latter form the Sheth, Mo & Tormen 
(2001) mass function is redshift independent. We probe to 
very high masses and very low abundances to illustrate the 
importance of high-accuracy merger kernels when consider- 
ing rare objects. The unmodified extended Press-Schechter 
merger kernel performs poorly, quickly resulting in a mass 
function shifted to low masses. The kernel with parame- 
ters identified by PCH08 performs much better, but also 
begins to underpredict the abundance of the most mass ha- 
los at high-redshift. The kernel using parameters obtained 
in this work by fitting conditional mass functions performs 
extremely well, remaining remarkably close to the expected 
Sheth, Mo & Tormen (2001) mass function out to Z = 4, al- 
though statistically significant differences can be detected. 
The kernel with parameters selected to best solve Smolu- 
chowski's equation (red line) also performs remarkably well. 
In particular, it is the most successful at match the evolution 
of the most massive halos in the mass function, as would be 
expected from Fig. 2. It performs somewhat worse than the 
kernel fit to conditional mass functions at low masses as also 
expected. 



tended Press-Schechter kernel). Finally, the red line shows 
results from the kernel which provides the best solution to 
Smoluchowski's equation. "Best" here is in the sense defined 
by eqn. (14), which tries to solve Smoluchowski's equation 
most accurately where the rate of change of the mass func- 
tion per halo is largest. Not surprisingly, therefore, this gives 
the most accurate solution for massive halos [z >1), but ac- 
tually gives a worse solution to Smoluchowski's equation at 
low masses than do kernels fit to conditional mass functions. 
This simply reflects the relative importance given to differ- 
ent mass ranges by each constraint and the fact that the 
particular functional form of the merger kernel chosen does 
not permit a precise solution to Smoluchowski's equation, 
thereby forcing a compromise solution to be adopted. 

Figure 3 (which has the same format as Figure 1 of 
PCH08) shows conditional mass functions from the Millen- 
nium Simulation (black histograms) for three different final 
halo mass ranges and for four different redshifts. Overplot- 
ted are conditional mass functions estimated from merger 
trees built using different merger kernels. The failure of the 
unmodified extended Press-Schechter kernel (green lines) is 
readily apparent, while blue lines — which use kernels con- 
strained by fitting to these same conditional mass functions 
in this work (solid lines) and PCH08 (dotted lines) — are 
vastly improved matches as expected. The red line indi- 



4 DISCUSSION 

We have described how Smoluchowski's equation, which gov- 
erns any mass-conserving binary coagulation process, may 
be used to provide constraints on halo merger rates. These 
constraints are complementary to those obtained by fitting 
to conditional mass functions from N-body simulations since 
they span a wide dynamical range of halo masses. Using 
these methods, and the modified extended Press-Schechter 
algorithm described by PCH08, we have identified merger 
kernels which best-fit conditional mass functions from the 
Millennium Simulation and which best solve Smoluchowski's 
equation. While the functional form of the PCH08 merger 
rate does not permit an exact solution of Smoluchowski's 
equation, we are able to find solutions which greatly improve 
the accuracy of merger trees. Using these best-fit solutions 
we are able to evolve the dark matter halo mass function 
from Z — to Z — 4 with remarkably high accuracy, par- 
ticularly for the most massive objects. 

Ideally, we would identify a functional form for the 
merger kernel which permits a precise solution to Smolu- 
chowski's equation, while simultaneously producing progen- 
itor masses functions consistent with the available N-body 
data. In fact, a perfect solution should agree with N-body 
data when compared using any statistic (e.g. distributions 
of most massive or second most massive progenitors). In re- 
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Figure 3. Comparison of progenitor mass functions measured from the Millennium Simulation (black histograms; taken from Cole et 
al. (2008)) and those constructed using merger trees for a variety of merger kernels, (Go, 71, 72). The green lines shows results for the 
standard extended Press-Schechter merger kernel, (1, 0, 0), while the red line shows results using the merger kernel which best solves the 
Smoluchowski equation. The solid blue curve is the best fit to the Millennium Simulation found in this work, while the dotted blue line 
uses the best fit from PCH08. For all calculations, we use the Sheth, Mo & Tormen (2001) mass function in Smoluchowski's equation. 
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Figure 4. Left-hand panel: The dark matter halo mass function shown at Z = 0, 0.5, 1, 2 and 4 (from right to left). The analytic Sheth, 
Mo & Tormen (2001) mass function is shown by black lines. Coloured lines show the result of evolving the mass function from Z = 
to higher redshifts using merger trees with different merger kernels, (Go, 71,72)- The green lines show results for the standard extended 
Press-Schechter merger kernel, (1,0,0), while the red lines show results using the merger kernel which best solves the Smoluchowski 
equation. The solid blue lines are for the the best fit to the Millennium Simulation conditional mass functions found in this work, while 
the dotted blue lines use the best fit to conditional mass functions from PCH08. 



ality, we do not know of such a function and, as mentioned 
in §1, only a handful of analytic solutions to Smoluchowski 's 
equation are known. In that case, the search for the "best" 
kernel requires making some decision about what are the 
most important statistics to fit and accepting that some de- 
gree of compromise in matching them is unavoidable. We 
believe that the most important statistics to match are the 
evolution of the overall mass function (which is extremely 
well constrained from N-body simulations) and progenitor 
mass functions. When considering the process of structure 
and galaxy formation it is these statistics which control, to 
first order, the number of galaxies able to form at a given 
point in cosmic history and how those galaxies were formed. 

Practically, these modified merger rates should prove 
extremely useful in constructing high-accuracy merger trees 
for use in studies of structure formation, reionization and 
galaxy formation. If even greater accuracy is required, a 
different functional form for the merger kernel must be 
adopted^ and the techniques described in this work ap- 
plied to constrain its parameters. As noted in §1 only a 
handful of analytic solutions to Smoluchowski's equation are 
known, none of which provide useful solutions for dark mat- 
ter halo merger rates in cold dark matter Universes. Never- 
theless, a sufficiently general parametrization of the merger 
kernel should permit arbitrarily accurate solutions to Smolu- 
chowski's equation. If used in a binary split merger tree algo- 
rithm such a kernel should allow for an arbitrarily accurate 



^ The functional form of PCH08's modification of the merger ker- 
nel can be considered to be the first-order terms in a Taylor expan- 
sion of In G — adding higher-order terms would be straightforward. 



evolution of the halo mass function (limited only by the ac- 
curacy of the merger tree construction algorithm) . 

This work does not provide any further physical in- 
sight into dark matter halo merger rates — at least not di- 
rectly. One might hope that physical insight into the "micro- 
physics" (to make an analogy with other areas in which 
Smoluchowski's equation is applied, e.g. polymer growth) of 
dark matter halo merging might point towards a functional 
form for the kernel. Until such insight is gained, the meth- 
ods described herein provide a practical method for rapid 
construction of high-accuracy merger trees. 
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APPENDIX A: RELATION BETWEEN MERGER TREE SPLIT PROBABILITIES AND 
SMOLUCHOWSKI EQUATION MERGER KERNEL 

Smoluchowski's equation is normally employed to take an existing distribution of masses and evolve it forwards in time subject 
to a specified merger kernel. Both creation and destruction of halos must be considered in this case. However, we can also apply 
Smoluchowski's equation to the reverse process, taking a distribution of halos of specified mass and evolving them backwards 
in time. In this case, we have a series of mass-conserving binary splits reminiscent of merger trees. The split probability is 
related to the creation term in Smoluchowski's equation. 

PCH08 describe a merger tree binary split algorithm which utilizes a modification of the extended Press-Schechter split 
probability distribution. Here, we wish to cast that same algorithm in the terms used in Smoluchowski's equation. 

Assuming a Sheth, Mo & Tormen (2001) mass function, Smoluchowski's equation is 

Usmt(z) = - J n SM T(z') n SMT(z - z')q(z', 2 - z')dz - j n S MT(z)nsMT(z')q(z,z')dz', (Al) 

where usmt(z) is the Sheth, Mo & Tormen (2001) mass function and hsmt(z) is the rate of change of that mass function. 

Applying Smoluchowski's equation in reverse, we find that in a binary split algorithm, the split probability distribution 
for a single halo of mass z to have a progenitor of mass z' (and, therefore, a second progenitor of mass z — 2') is given by 

,. = n SMT {z')ns^{z-z') q{ ^ z _, y 

The PCH08 split probability function is that derived from the Press-Schechter mass function and extended Press-Schechter 
techniques: 

. , . nps(z')np S {z - z) , , , 

P ep nt(z ; z) = K —t — -A- -<7 c ps(z ,z-z), (A3) 

nps(z) 

where q' cPS (zi,Z{ — z\) = g e ps(«i,«f — z{)G{o\/oi , Si/at) and G(ai/(Tf,5f/af) is PCH08's multiplicative modifier of merging 
rates. 

From this, we can deduce that the rate at which systems of mass 21 and 22 merge together is just 

R(zi, 22) = n S MT(2i + 2 2 )P S piit(2i; 21 + 22) = nsMT ( Zl + 2 ^ n P s(2i)wps(22)gcPs(zi, 22) (A4) 

nps(2i + 22 j 

Therefore, 

ysMT(z) = p-r- / nps(2 )n PS (2 - 2 )<7ops(2 ,2 - 2 )dz - / — ■ — — n PS (2)nps(2 )g oPS (2, 2 )d2 , (A5) 

nps(z) 2 J Q J () n PS (2 + 2') 

or, in terms of the rate of change per halo, 
?/smt(2) _ 1 



-/ n PS (2 )nps(z- 2 )g cPS (2 ,2-2 )dz - I — — ■ — -np S {z)np S (z )q e p S (z, z )dz 

2 Jo J nsMT{z)nps{z + 2' 



usmt{z) nps(z) 

This is the form of Smoluchowski's equation that we utilize throughout this work. 



■(A6) 
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